Magnetic susceptibility of QCD at zero and at finite temperature from the lattice 
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The response of the QCD vacuum to a constant external (electro)magnetic field is studied through 
the tensor polarization of the chiral condensate and the magnetic susceptibility at zero and at finite 
temperature. We determine these quantities using lattice configurations generated with the tree- 
level Symanzik improved gauge action and Nf = 1 + 1 + 1 flavors of stout smeared staggered 
quarks with physical masses. We carry out the renormalization of the observables under study 
and perform the continuum limit both at T > and at T = 0, using different lattice spacings. 
Finite size effects are studied by using various spatial lattice volumes. The magnetic susceptibilities 
Xf reveal a diamagnetic behavior; we obtain at zero temperature \ u — —(2.08 ± 0.08) GeV -2 , 
Xd = -(2.02 ± 0.09) GeV~ 2 and X s = -(3.4 ± 1.4) GeV~ 2 for the up, down and strange quarks, 
respectively, in the MS scheme at a renormalization scale of 2 GeV. We also find the polarization 
to change smoothly with the temperature in the confinement phase and then to drastically reduce 
around the transition region. 



INTRODUCTION 



An external (electro)magnetic field is an excellent 
probe of the dynamics of the QCD vacuum. Strong 
magnetic fields affect fundamental properties of QCD 
like chiral symmetry breaking and restoration, deconfinc- 
ment, the hadron spectrum or the phase diagram, just 
to name a few. Chiral symmetry breaking has long been 
known to be enhanced by magnetic fields at zero temper- 
ature, signalled by an increasing chiral condensate (see 
e.g. Ref. [1]). The particle spectrum may undergo dras- 
tic changes (see e.g. the ongoing discussion in Refs. [2- 
4]) with some strong decay channels becoming unavail- 
able and others opening up. The transitions at non- 
vanishing temperature related to chiral symmetry break- 
ing and deconfinement are also affected by the magnetic 
field B. The phase diagram of QCD in the temperature- 
magnetic field plane was determined recently in lattice 
simulations [5-7] by analyzing the dependence of the chi- 
ral condensate and of other observables on B, with the 
main result that the transition temperature T c decreases 
with growing B 1 and the transition remains an analytic 
crossover just as at B = [10]. These effects are relevant 
in several physical situations as strong magnetic fields 
are expected to play a significant role, e.g., in early cos- 
mology [11], in non-central heavy ion collisions [12] and 
in dense neutron stars [13]. 



* Corresponding author, gergely.endrodi@physik.uni-r.de 
1 Employing physical quark masses in the simulation and extrapo- 
lating the results to the continuum limit, as was done in Refs. [5— 
7], proved to be essential. Studies where these ingredients are 
missing produce qualitatively different results, namely an in- 
creasing T C (B) function [8, 9]. A possible explanation for this 
discrepancy and a comparison to effective theories was given re- 
cently in Ref. [7] . 



Another fundamental characteristic of the QCD vac- 
uum is the response of the free energy density (which at 
zero temperature is the vacuum energy density) to mag- 
netic fields, 
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where Z is the partition function of the system and V 
the (three-dimensional) volume. Due to rotational invari- 
ance the independence of / is to leading order quadratic, 
characterized by the magnetic susceptibility of the QCD 
vacuum. 
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which is a dimensionless quantity (here e > denotes the 
elementary charge). A positive susceptibility indicates a 
decrease in / due to the magnetic field, that is to say, 
a paramagnetic response. On the other hand £ < is 
referred to as diamagnetism [14]. Clearly, the sign of £ is 
a fundamental property of the QCD vacuum. 

In the functional integral formalism of QCD the sus- 
ceptibility is readily split into spin- and orbital angular 
momentum-related terms, according to 
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with contributions from each quark flavor / with electric 
charge qf and mass nif. For a constant magnetic field 
B = F xy in the positive z direction, 
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£j is given by an analogous expression with <j xy replaced 
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by a generalized angular momentum also present for spin- 
less particles, cf. Eq. (A7) of Appendix A. Eq. (4) consti- 
tutes an important relation which, to our knowledge, has 
not been recognized previously in this context. Its deriva- 
tion from the quark determinant and the corresponding 
Dirac operator is given in Appendix A. 

In the present paper we concentrate on the spin con- 
tributions, and thus the expectation value of the tensor 
polarization operator ipfa^ipf- To leading order this is 
proportional to the field strength and thus can be written 
as [15] 

{$f(Tx V i)f) = QfB ■ (4>f*l>f) ■ Xf = QfB ■ T f , (5) 

where Oi/jfipf) is the expectation value of the quark con- 
densate (see Eq. (9) below for its definition). Corrections 
to the right hand side are expected to be of 0(B 3 ), so 
that Lorentz invariance is maintained. In the literature 
X/ is referred to as the magnetic susceptibility of the con- 
densate (for the quark flavor /). In what follows we will 
also use the term "magnetic susceptibility". Again we 
stress that it constitutes only one of the two contribu- 
tions to the total susceptibility. We also define the ten- 
sor coefficient Tf as the product of the condensate and 
the magnetic susceptibility. Both quantities will depend 
on the temperature T at which the expectation values 
of Eq. (5) are determined. At finite quark masses it is 
advantageous to work with r/ instead of Xf f° r reasons 
related to renormalization (see below). 

The magnetic susceptibility X/, in the context of QCD, 
was first introduced in Ref. [15]. Since then its experi- 
mental relevance has been growing steadily. In particu- 
lar, this quantity appears in the description of radiative 
D s meson transitions [16], of the anomalous magnetic 
moment of the muon [17] and of chiral-odd photon dis- 
tribution amplitudes [18, 19]. Moreover, vector-tensor 
two-point functions at zero momentum are related to the 
magnetic susceptibility [20]. 

Since Xf ac ts as an input parameter in various strong 
interaction processes [21], a high-precision determination 
of its value is of importance. In the past, the mag- 
netic susceptibility has been calculated using QCD sum 
rules [22-24], in the holographic approach [25, 26], us- 
ing the operator product expansion [27], in the instan- 
ton liquid model [28] and in low-energy models of QCD 
like the quark-meson model and the Nambu-Jona-Lasinio 
(NJL) model [29]. The numerical value of Xf was a ls° 
determined recently on the lattice in the quenched ap- 
proximation of two- [30] and of three-color QCD [31], in 
both cases without renormalization. We mention that 
the quenched approach can lead to large systematic er- 
rors at strong magnetic fields [32]. 

In this paper we determine Tf(T) and X/(T) for a wide 
range of temperatures around the transition region be- 
tween the hadronic and the quark-gluon plasma phases 
and at T = 0. We apply fully dynamical lattice simu- 
lations, i.e. both the fermionic degrees of freedom and 
the external field are taken into account in the gener- 



ation of the gauge ensembles. We perform the renor- 
malization of the tensor coefficient and carry out the 
continuum extrapolation using results obtained at dif- 
ferent lattice spacings. One main result will be that the 
tensor coefficients at T = are negative, indicating the 
spin-diamagnetic nature of the QCD vacuum. Moreover 
we observe that Tf decreases around the QCD crossover 
temperature similarly to other order parameters like the 
condensate. 

This paper is organized as follows. We define the lat- 
tice implementation of the magnetic field and the observ- 
ables in Sec. 2 and discuss their renormalization in Sec. 3. 
The multiplicative renormalization is carried out pertur- 
batively; the determination of renormalization constants 
is detailed in Sec. 4. After a brief summary of the simula- 
tion setup in Sec. 5 we present the results in Sec. 6 for the 
tensor coefficients and in Sec. 7 for the susceptibilities, 
before we conclude. 



2. MAGNETIC FIELD AND OBSERVABLES 

We study the effect of an external magnetic field 
B on the expectation value of the tensor polarization, 
Eq. (5). To realize such an external field on the lattice 
we implement the continuum U(l) gauge field satis- 
fying d x A y — d y A x = B using space-dependent complex 
phases [5, 8, 33, 34] in the following way, 

Uy (n) = e " l2q f Bn *, 
u x (N x - l,n y ,n z ,n t ) = e -- 2 «/^n K } 

u x (n) = l, n x ^N x -l, 
U v [n) = l, v^{x,y}, 

where the sites are labeled by integers n = 
(n x , n yi n z , n t ), with n v = . . . N v — 1 and a is the lat- 
tice spacing. This prescription for the links corresponds 
to a covariant derivative for the flavor / of the form 
D^j = + iqfA^ + iA^T a 2 . This discretization satis- 
fies periodic boundary conditions in the spatial directions 
and ensures that the magnetic flux across the x — y plane 
is constant. It is well known that the magnetic flux in a 
finite volume is quantized [35, 36], which on the lattice 
implies 

q B-a 2 = ^f, N b eZ, 0<N b <N x N v , (7) 

where q is the smallest charge in the system, in our case 
q = q c i = q s = — e/3. Due to the periodicity of the 
links of Eq. (6) in Nb with period N x N y , one expects 



2 Note that we do not include in the action the corresponding 
photon kinetic term F^F^ j '4 = B 2 /2. This means that in the 
discussion we will never encounter B alone but only the combi- 
nation QfB ~ eB. 
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lattice artefacts to become large if iV& > N x N y /A. In the 
following we use lattices with N x = N y — N z = N s . 

We consider three quark flavors u, d and s. Since the 
charges and masses of the quarks differ we have to treat 
each flavor separately; q u = —2q c i = —2q s . We assume 
m u = rrid m s . The partition function in the staggered 
formulation then reads, 

Z = J VUe- f3S * Y[ [det M(U, q f B, m/)] 1/4 , (8) 

with M(U,qB,m) = lp(U,qB) + ml being the fermion 
matrix and j3 = 6/g 2 the gauge coupling. The exact form 
of the action we use is described in Refs. [37, 38], and 
further details of the simulation setup are given in Sec. 5. 
Since the external field couples directly only to quarks, B 
just enters the fermion determinants, through the U(l) 
links of Eq. (6). The volume of the system is given as 
V = (aN s ) 3 and the temperature as T = (a7V t ) _1 . 

In this formulation the expectation value of the quark 
condensate for the flavor / can be written as 

55 V-^T = & {*M-\U, qf B,m f) ) . (9) 

Likewise, the expectation value of the tensor Dirac struc- 
ture reads, 

(fav^f) = ^ (trM-^l/.g/B.m/)^) . (10) 

The staggered realization of this operator is detailed in 
Appendix C. 

At this point a few comments regarding the sign of 
the expectation values in Eq. (5) are in place. In con- 
tinuum calculations a negative sign for the condensate 
is customary, see e.g. Ref. [24], in contrast to our con- 
vention in Eq. (9). This sign convention applies for any 
fermionic bilinear expectation value, therefore it does not 
affect the sign of Xf, but only that of r/. Further pos- 
sible differences in the sign can arise from the definition 
of a^y and from that of the U(l) part of the covariant 
derivative. We note that our notation is consistent with 
that of Ref. [24] in terms of cx M „, but differs by a minus 
sign in the covariant derivative (see the paragraph below 
Eq. (6)), implying an overall relative minus sign of Xf- 



3. RENORMALIZATION 

In order to determine the continuum limit of the ob- 
servables defined in Eqs. (9) and (10), their renormaliza- 
tion has to be performed. The quark condensate (at finite 
mass) is subject to additive and multiplicative renormal- 
ization, due to the divergent terms in the free energy den- 
sity / of Eq. (1) and in the bare mass m/. The former 
divergence is (to leading order) quadratic in the cutoff 
1/ot [39]. Therefore, the bare observable can be written 



<^/) (B,T) = ~ (V^/) r (B,T) + Csm f /a 2 + ..., 

where Zg is the renormalization constant of the scalar op- 
erator and the ellipses denote subleading (logarithmic) 
divergences in a. Here the superscript r indicates the 
renormalized observable. The divergences in (/tpy) de- 
pend neither on the temperature nor on the external 
field 3 . Therefore, in mass-independent renormalization 
schemes, C,s and Z$ are just functions of the gauge cou- 
pling. The conventional way to cancel the additive diver- 
gences is to consider the difference, for example, between 
the condensate at T and at T = 0. 

The situation is somewhat different for the tensor po- 
larization. As a calculation in the free theory shows, 
an additive divergence of the form qfBrrif log(m 2 a 2 ) ap- 
pears in NxTuvip) (see Appendix B). This divergence van- 
ishes in the chiral limit (or at zero external field) and is 
not related to the multiplicative divergence of the tensor 
operator to which we will return below. Altogether the 
bare observable can thus be written as 

<^/tvV/>(B,T) 
= ^- (^f^xy^fY (B, T) + ( T q f Bm f log(m 2 f a 2 ) + 

(12) 

where Zt is the renormalization constant of the ten- 
sor operator (its perturbative determination is detailed 
in Sec. 4) and £t the coefficient of the divergent loga- 
rithm. Both are independent of T and B (and in mass- 
independent schemes of m/). In Eq. (12) the ellipses 
denote finite terms. In the free theory we calculate 
Ct(9 = 0) = 3/(4tt 2 ) (see Appendix B). 

From these considerations it is clear that the mag- 
netic susceptibility Xfi being proportional to the ratio 
of Eq. (12) over Eq. (11), at non-vanishing quark mass 
contains additive divergences which depend both on T 
and on B (and also on the quark flavor /). This means 
that these singular contributions cannot be removed by 
subtracting the same operator, measured at different T 
or B (or flavor /). 

Therefore, in the following we consider the tensor coef- 
ficient Tf defined in Eq. (5). We notice that the operator 
1 — rrifd/dmf eliminates the logarithmic divergence and 
thus can be used to define an observable with a finite 
continuum limit, 

r; s (l - m,^j r r Z T ^ r f Z T - . (13) 

At finite quark mass this is one possible prescription to 



3 For a detailed argumentation about the absence of B-dependent 
divergences in the condensate see Ref. [5] and references therein. 
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cancel the additive logarithmic term. It has the advan- 
tages that the chiral limit of t/ is left unaffected, and 
that, together with the logarithmic divergence, scheme- 
dependent finite terms also cancel in this difference (see 
Eq. (B7)), such that the scheme- and renormalization 
scale-dependence of t/ resides solely in Z^. 

Since the subtracted divergence is independent of the 
temperature, we are able to determine t^ 1v at zero tem- 
perature where we systematically study the dependence 
of Tf on m/ and a, and then to perform the subtraction 
at nonzero temperatures as well. As we will see, the sub- 
traction in Eq. (13) amounts to a 5 — 10 per cent effect 
for the lattice spacings we use. 

4. PERTURB ATIVE EVALUATION 

4.1. Fermion propagator and fermion bilinears 

The discretization of the continuum bilinear operators 
in the staggered formulation of lattice QCD is detailed 
in Appendix C. Using this discretization we compute the 
one-loop correction to the fermion propagator in order to 
obtain its renormalization constant, an essential ingredi- 
ent for the renormalization of fermion operators defined 
in Eqs. (C7) and (C8). The one-loop Feynman diagrams 
that enter our two-point amputated Green function cal- 
culation for the propagator are illustrated in Fig. 1. For 
the algebraic operations involved in evaluating Feynman 
diagrams, we make use of our symbolic package in Math- 



ematica. The required procedure for the computation of 
a Feynman diagram can be found in Ref . [40] . 




1 2 

FIG. 1: One-loop diagrams contributing to the fermion prop- 
agator. Wavy (solid) lines represent gluons (fermions). 

We provide the total expression for the inverse fermion 
propagator S~ 1 (p), computed up to one-loop in pertur- 
bation theory. Here we should point out that, for dimen- 
sional reasons, there is a global prefactor 1/a multiplying 
our expressions for the inverse propagator, and thus, the 
O(a ) correction is achieved by considering all terms up 
to 0(a 1 ). We have computed S~ 1 (p) for general values 
of the gauge parameter a, the stout smearing parameters 
lui,uj2 (since we perform two stout steps, we have kept 
these different, so that our results are also applicable for 
the single stout smearing case), the bare quark masses 
and the external momentum. We have obtained results 
using 10 sets of values for the Symanzik improvement co- 
efficients of the gluon action (shown in Ref. [40]). In pre- 
senting our result for S 1-1 ^) below, we use the Landau 
gauge for conciseness; the quantities Cj are numerical co- 
efficients that depend on the gluon action we choose and 
on the stout smearing parameters. 
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N c = 3 is the number of colors and q\ and q 2 are the 
external momenta. Next we study the one-loop per- 
turbative Green functions of the tensor fermionic opera- 
tor. The one-particle irreducible Feynman diagrams that 
enter the calculation of the above Green functions are 
shown in Fig. 2, and include up to two-gluon vertices ex- 
tracted from the operator (the cross in the diagrams). 
The appearance of gluon lines attached to the operator 
stems from the product of gauge links in the operator 
definition (see Eq. (C4)). 




X 



FIG. 2: One-loop diagrams contributing to the bilinear op- 
erators. A wavy (solid) line represents gluons (fermions). A 
cross denotes the insertion of a Dirac matrix. 



4.2. Quark field and quark bilinear renormalization 
constants in the RI'-MOM scheme 



The renormalization constant (RC) can be thought of 
as the link between the matrix element, regularized on 
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the lattice (without upper index), and the renormalized 
continuum counterpart (with upper index r). For the 
fermion field, the scalar and the tensor operator this is 
written as 



ZoO, 



or. 



Zrr On 



(15) 



The RCs of lattice operators are necessary ingredients 
in the prediction of physical probability amplitudes from 
lattice matrix elements. In this section we present the 
multiplicative RCs, in the RI'-MOM scheme, of the quark 
field (Z q ) and the quark bilinear operators (Z$ and Z^). 

In the RI'-MOM renormalization scheme the forward 
amputated Green function A(p), computed in the chiral 
limit and at a given (large Euclidean) scale p 2 = /i 2 , 
is imposed to be equal to its tree- level value. The RCs 
are computed at arbitrary values of the renormalization 
scale. 



^i-ioopO 3 

Al-loop(P 



/i, m 



= 0) = S u .l c (n,m = 0)Z q (fi), 



fi, m = 0) = A t s rcc (^, m = Q)^ q (/i) Z s l (fi), 



where s = S corresponds to the scalar and s = T to 
the tensor bilinear operator. Moreover, here /j, is the 
renormalization scale, S^ c is the tree- level result for the 
inverse propagator, and Af rcc is the tree-level value of the 
Green function for the operator under study. 



The expressions we obtain using our one-loop results 
discussed in the previous subsection are shown here only 
for the tree-level improved Symanzik gauge action. We 
allow for a general gauge characterized by the parame- 
ter a (Landau gauge: a = 0, Feynman gauge: a = 1) 
and stout smearing parameters. As discussed above, we 
have employed different parameters for the two smear- 
ing steps; in fact, we have also kept the parameters of 
the action's smearing procedure (w^, uja 2 ) distinct from 
the parameters of the operator smearing (woj,wo 2 ). In 
all expressions the systematic error (coming from an ex- 
trapolation to infinite lattice size of our numerical loop- 
integrals) is smaller than the last digit we present. The 
results read: 
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4.3. Conversion to the MS scheme 

Here we provide the expressions for the RCs in the 
MS continuum scheme, using conversion factors adapted 
from Refs. [41, 42]. These conversion factors do not de- 
pend on the regularization scheme (and, thus, they are 
independent of the lattice discretization) , when expanded 
in terms of the renormalized coupling constant. However, 
expressing them in terms of the bare coupling constant 
in general introduces a dependence on the action. To 
one-loop order, the renormalized and the bare couplings 
are nevertheless equal, leading to 
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(20) 

We note that the RCs from the lattice scheme to the 
MS scheme arc independent of the gauge parameter a. 
as they should be. 

The above conversion factors refer to the NDR (Naive 
Dimensional Regularization) version of MS (see, e.g., 
Ref. [43]). Other possible modified minimal subtrac- 
tions are related to NDR via additional finite renormal- 
izations. In particular, to compare our result for Zt 
with the one in Ref. [44], which is given in the "DREZ" 
scheme (for the Wilson gauge action and without stout 
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smearing), we must divide our NDR result by a factor 
(1 + Ci?<7 2 /(167r 2 )). The two results are in perfect agree- 
ment. 



5. SIMULATION SETUP 

For our measurements we used the gauge ensembles 
of Refs. [5, 7] augmented by additional new ensembles. 
All configurations were generated with the tree-level im- 
proved Symanzik gauge action and stout smeared stag- 
gered fermions, at physical quark masses. We use lat- 
tices at both T = and at T > 0, at various values 
of the external magnetic field. We employ two steps of 
stout smearing with parameter u>a — ujo — 0.15 both in 
the action and in the operators. The zero temperature 
ensembles consist of 24 3 x 32, 32 3 x 48 and 40 3 x 48 lat- 
tices at five different lattice spacings, while at finite tem- 
perature we carried out measurements on lattices with 
N t — 6, 8 and 10, to enable a continuum limit extrap- 
olation. We studied finite volume effects on Nt = 6 
lattices, using three different aspect ratios. The light 
(m u = m d = m u d) and strange (m s ) quark masses are set 
to their physical values, along the line of constant physics 
(LCP) as m ud = m ud (f5), and m s /m ud = 28.15. The 
LCP was determined by keeping fx/M n and Jk/^k 
physical, and the lattice scale is set using More de- 
tails about the lattice action, the determination of the 
scale and the LCP, and the lattice ensembles can be found 
in Refs. [5, 37, 38]. At each temperature and external 
magnetic field we measured the observables of interest 
on 0(100) thermalized configurations which were sepa- 
rated by 5 trajectories to reduce autocorrelations. The 
measurements were carried out using the noisy estimator 
method, with 20-40 random vectors. 

For our choice of the smearing parameters we obtain 
for the scalar and tensor renormalization constants, 
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from Eqs. (17), (18) and (20). We define the coupling g 
in the "E" scheme [45], using the plaquette expectation 
value, 



9% 
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(trUb) 



(22) 



which is found to be 10 — 20% larger than the bare cou- 
pling g 2 . We compute c perturbatively for the tree- 
level improved Symanzik gauge action and obtain c = 
0.183131340(2) • C F , thereby confirming Ref. [46]. 

We allow for a systematic error of 50% in 1 — Z^ s 
for the effect of higher order terms in the perturbative 
calculation. 



6. RESULTS 

We measure the tensor polarizations as functions of 
the external field at various temperatures for the three 
different flavors. We observe that (ti> u (X xy' , Pu) is negative, 
indicating that Xu < 0, in accordance with Ref. [24] and 
the discussion about the sign convention below Eq. (10). 
Whether this corresponds to a para- or a diamagnetic 
response will be discussed in Sec. 7. 
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FIG. 3: Minus the bare tensor polarization (upper panel) and 
the bare condensate (lower panel) for the up quark, for three 
temperatures on the Nt — 6 lattices. 

In the upper panel of Fig. 3 we show the bare tensor po- 
larization as a function of the magnetic field for Nt = 6. 
We confirm the linear trend to leading order in B, in 
agreement with Ref. [30]. However, the slope at small B 
is also observed to change significantly with temperature. 
We find that nonlinear effects are always below 5% for 
magnetic fields eB < 0.2 GeV 2 and they reduce as the 
temperature decreases. In the lower panel of Fig. 3 we 
also show how the bare condensate itself changes with B 
for different temperatures. We observe that the depen- 
dence of the condensate on B varies strongly with the 
temperature in the transition region. This behavior was 
found to be the reason for the decrease of the chiral tran- 
sition temperature with growing B, and was investigated 
in detail in Refs. [5-7]. We study finite volume effects 
at one temperature T = 141 MeV, for N t = 6 ensembles 
with N s — 16, 24 and 32, see Fig. 3. The largest lattice 
corresponds to a linear extent of 7 fm. Since we see no 
deviation for the tensor polarization or the condensate 
between the different volumes, we conclude that finite 
size effects are smaller than our statistical errors. 
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Next, we concentrate on the leading linear trend in 
{^f^xyfpf), i- e - on the slope characterized by the tensor 
coefficient Tf, as defined in Eq. (5). We perform the mul- 
tiplicative renormalization of Tf according to Eq. (13), 
using the tensor renormalization constant, Eq. (21) to the 
MS scheme at a renormalization scale /i = 2 GeV. The 
dependence of the results on the renormalization scale fi 
is found to be mild, as can be seen below. 

We measure Zt ■ Tf at zero temperature for several lat- 
tice spacings and quark masses. Here we fix the strange 
quark mass to its physical value and tune only the light 
mass such that R = m uf j/™^ ys varies between 0.5 and 
28.15. For the latter ratio all three quarks have equal 
masses. (Note that these measurements are also fully dy- 
namical and no partial quenching is applied.) In Fig. 4 
we plot minus the tensor coefficient for the up quark as 
a function of R for five different lattice spacings. Moti- 
vated by the behavior of the tensor coefficient in the free 
case, Eq. (B6), and by the scaling properties of the action 
we use, we consider the following fit function for Z^Tf. 

c f0 + c fl R + c f2 Rlog(R 2 a 2 ), c fi = cf) + c^a 2 . (23) 

Here a is to be understood in units of GeV -1 . This form 
describes the data very well; we obtain x 2 /dof. < 1.5 for 

f i) 

both the up and down flavors. The fitted values for c\ ? 
are listed in Table I. We remark that the coefficients of 
the logarithms, c^/m* 8 = 0.055(5) and c$ /m^ s = 
0.072(6) are quite close to the free-field value of 3/(47r 2 ) 
(see Appendix B). We perform the fit both for all lattice 
spacings and for only the finest four lattices. Moreover 
we consider the inclusion of an R 2 term in the fit and 
vary the fit range to exclude points with largest masses. 
The difference between these fits is used to estimate the 
systematic error of this combined extrapolation. 




m „ / mP^ s 

ud ' ud 

FIG. 4: Mass dependence of the combination — Zt ■ t u in the 
MS scheme at renormalization scale /i = 2 GeV. The coef- 
ficient of the logarithmic divergence is determined by fitting 
the data by a lattice spacing-dependent function (solid lines). 

At zero quark mass the additive divergence is absent 



/ 


c (0) 


c (1) 


c (0) 


c (1) 

Cfl 


c (0) 


c (1) 


u 


-40.3 


3.8 


-2.1 


0.5 


0.19 


-0.03 


d 


-38.9 


2.8 


-2.5 


0.7 


0.25 


-0.07 



TABLE I: Central values for the fit parameters of Eq. (23) in 
units of MeV. 

and therefore, applying the combined fit, the continuum 
limit of the chiral limit of Z^Tf can be extracted (it equals 

the cJq parameter). This corresponds to the black point 
in Fig. 4. However, since we are interested in the ten- 
sor coefficient at physical quark masses, we now follow 
the scheme of Eq. (13), subtracting the logarithmic di- 
vergence. We apply the operator 1 — rrifd mj — 1 — RBr. 
which acting on the fit function of Eq. (23) yields 

t) =c f0 -2c f2 R. (24) 

As already emphasized in Sec. 3, the subtraction of the 
divergent term Tf lv does not affect the chiral continuum 
limit since it vanishes at mt = 0. Moreover, this sub- 
traction eliminates the scheme-dependent finite terms (cf. 
Eq. (B7)), making the conversion to the MS scheme triv- 
ial. 

For the strange quark we do not perform a similar anal- 
ysis with modified strange quark masses, but subtract 
the logarithmic divergence by using the fit parameters 
for the down quark and R = 28.15. We find that the de- 
pendence of the strange quark tensor polarization on the 
light quark masses is below a few per cent (1% for the 
coarsest and 4% for the finest lattice). Therefore this ap- 
proximation introduces errors smaller than those already 
present due to statistics and renormalization. 

After the subtraction, the renormalized tensor coeffi- 
cient rj has a well defined continuum limit even for fi- 
nite quark masses. We find that for physical light quark 
masses \t£™\ < 2.5 MeV for our range of lattice spac- 
ings. For the strange quark the divergent contribution 
is larger in magnitude, giving rise to larger errors due to 
this subtraction. 

Our final results for the zero temperature renormalized 
tensor coefficients in the MS scheme at a renormalization 
scale ix = 2 GeV are summarized in Table II. For the light 
flavors this Table contains the results both for physical 
quark masses and for the chiral limit. These values may 
be compared to the unrenormalized quenched SU(2) lat- 
tice result —r u d — 46(3) MeV of Ref. [30] and to a similar 
study in the quenched SU(3) theory, 52 MeV [31]. Our re- 
sults are in reasonable agreement with the QCD sum rule 
result 50(15) MeV of Ref. [24], which was calculated at 
/i = 1 GeV (note that the scale dependence of Tf is small 
due to its small anomalous dimension, see Eq. (21)). We 
also compare our results to the NJL and quark-meson 
model predictions of 69 MeV and 65 MeV [29], respec- 
tively, which were obtained at an even lower renormal- 
ization scale of ~ 0.6 GeV. We remark that the same 
authors obtain a lower value of 44 MeV in the renormal- 
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ized version of the quark meson model [29] . 



J 


TfL 


f 


error 


stat. 


mult. 


cont. 


scale 


total 


u 


phys. 


-40.7 


0.2 


0.3 


1.0 


0.8 


1.3 


chir. 


-40.3 


0.2 


0.3 


1.1 


0.8 


1.4 


d 


phys. 


-39.4 


0.3 


0.3 


1.1 


0.8 


1.4 


chir. 


-38.9 


0.3 


0.3 


1.3 


0.8 


1.5 


s 


phys. 


-53.0 


0.5 


0.8 


7.1 


1.1 


7.2 



TABLE II: Results and error budget for the renormalized ten- 
sor coefficients for physical quark masses (phys.) and in the 
chiral limit (chir.). Given are (in units of MeV) the errors 
related to statistics, the multiplicative renormalization, the 
combined continuum fit, the lattice scale and, finally, the to- 
tal error. 

Next, we use the fact that the t^ iv contribution is in- 
dependent of T to perform the additive renormalization 
of the tensor coefficient at finite temperatures. In the 
upper panel of Fig. 5 we plot — t£ as a function of the 
temperature for three lattice spacings. We perform a si- 
multaneous fit of the results for different lattice spacings 
to an ^-dependent spline function. This dependence is 
of the form Nf 2 , once again to reflect the scaling proper- 
ties of our lattice action. We can read off the continuum 
extrapolation at N^ 2 = 0, which is shown in the figure 
by the hatched yellow band. The systematic error of the 
continuum extrapolation is estimated to be 1 MeV based 
on our experience at T — (see Table II) and is added to 
the statistical error in quadrature. Moreover, the uncer- 
tainty in the determination of the lattice scale (for details 
see Ref. [38]) propagates into this result and gives rise to 
an additional systematic error of 2%. Since this latter 
error is uniform and does not influence the shape of the 
t t AT) curve, it is not included in the plot. 





i i i I i i 

- 


i I i i i I i 


i i I i i i I i i i 
continuum 


30 






— N,= 10 J 

— N,=8 








— n,=6 : 


20 








10 


7 MS, fj,= 2 

phys 


GeV 







— 1 1 1 1 1 1 


1 1 I I I 1 I 






120 


140 160 


180 200 



T (MeV) 

FIG. 5: Minus the renormalized tensor coefficient tJ(T) in 
the MS scheme at a renormalization scale /i — 2 GeV for 
three lattice spacings and the continuum extrapolation. 

In the same manner we determine the tensor coefficient 



for the down quark at T > 0, and obtain results which 
are within errors consistent with t£, just as was observed 
at T = 0. For the strange quark this procedure leads to 
a qualitatively similar temperature-dependence too. The 
dependence of d on the temperature in the transition 
region can be used to define a transition temperature at 
5 = 0. We determine the inflection point of r^ d (T) 
and obtain T c — 162(3) (3) MeV in the continuum limit. 
Here the first error combines the statistical error and the 
error coming from the continuum extrapolation, and the 
second one is due to the uncertainty in the lattice scale. 
In conclusion, the tensor coefficient acts as a quasi-order 
parameter for the chiral transition, and gives a similar 
transition temperature as the chiral condensate at B = 0, 
T c = 159(3) (3) MeV, cf. Refs. [5, 47]. 

Finally, we study the dependence of on the renor- 
malization scale fj, at T = 0. We carry out the analy- 
sis for a range of renormalization scales in the window 
1 GeV < /i < 4 GeV. We find a very mild dependence 
on [i such that the tensor coefficients remain within the 
total errors given in Table II. 



7. MAGNETIC SUSCEPTIBILITY 

We can translate the result for rj to the magnetic sus- 
ceptibility X/ of Eq. (5) using the (scale- and scheme- 
dependent) value of the quark condensate. We recall the 
Gell-Mann-Oakes-Renner relation, 

M 2 F 2 = {m u + m d ) • + . . . , (25) 

which, at zero external field and in the chiral limit, re- 
lates the light condensate I — u, d to the quark masses 
and to the pion mass and decay constant, with F — 
86.2(5) MeV [48]. We make use of a recent lattice de- 
termination [49] of the quark masses in the MS scheme 
at_ fi = 2 GeV, m u + m d = 6.94(13) MeV, to extract 
(tpitpi) = (269(2) MeV) 3 . (We mention that multiplying 
the lattice bare mass along the LCP [38] and the inverse 
of the scalar renormalization constant of Eq. (21), we get 
a compatible value for the renormalized quark mass in 
the MS scheme, albeit with large uncertainties.) For the 
strange condensate we employ the QCD sum rule predic- 
tion [50], (ip s i/} s ) I (i'l'ipi) — 0.8(3). Using these values 
for the quark condensates, the zero-temperature mag- 
netic susceptibilities at physical quark masses are calcu- 
lated as 

Xu = -(2.08 ±0.08) GeV~ 2 , 
MS, [i = 2 GeV : Xd = -(2.02 ± 0.09) GeV" 2 , (26) 
Xs - -(3.4 ±1.4) GeV~ 2 . 

The magnetic susceptibilities at different values of [i can 
be obtained by running down with the ratio of renormal- 
ization constants Z^ s /Z^- s . Using the four-loop running 
to /i = 1 GeV one has to multiply the above values by 
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r = 1.49(7). We remark furthermore that running down 
with Zcf s to a renormalization scale of n = 1 GeV we 
obtain (rfiipt) = (245(5) MeV) 3 . 

Our results in Eq. (26) are in good agreement with the 
QCD sum rule calculations 4 summarized and updated in 
Ref. [24]: X i = -2.11(23) GeV~ 2 at fi = 2 GeV, and also 
compare well with the vector dominance estimate of xi = 
—2/m 2 p s» —3.3 GeV~ 2 . We remark that for the strange 
susceptibility, QCD sum rules predict \s ~ Xi [51], which 
is somewhat smaller than our result in Eq. (26). 

Comparing the temperature-dependence of the light 
tensor coefficient (Fig. 5) and that of the light quark 
condensate from Ref. [7], we conclude that the ratio of 
the two renormalized observables is compatible with a 
constant, resulting in a magnetic susceptibility xi(T) de- 
pending only weakly on the temperature, at least for tem- 
peratures T < 170 MeV. Moreover, we remark that since 
Xf is given in terms of the chiral condensate (which has a 
large anomalous dimension), the magnetic susceptibility 
has a stronger scale dependence than rj. 

As anticipated in the introduction, the magnetic sus- 
ceptibility Xf of the condensate is intimately connected 
to the spin contribution £ s to the total magnetic sus- 
ceptibility. Using this equivalence (which we prove in 
Appendix A), one sees that with our sign conventions 
Xf > corresponds to paramagnetism and Xf < to 
diamagnctism. Thus we conclude that the response of 
the QCD quark condensate to external magnetic fields is 
in its nature diamagnetic. 



8. CONCLUSIONS 

In this paper we studied the response of the QCD vac- 
uum to a constant external magnetic field at zero and at 
finite temperature. We determined the tensor polariza- 
tions of the quark condensates for various temperatures 
and external fields. We observed that the polarization of 
the flavor / at a temperature T is a linear function of B 
for fields eB < 0.2 GeV 2 , with a coefficient t/(T), defined 
in Eq. (5) . The renormalization of this tensor coefficient 
requires two steps. The additive divergences (which are 
present for finite quark masses) were fitted explicitly at 
T = and then subtracted using the operator 1 — md m , 
at T = and at T > 0. The multiplicative renormaliza- 
tion was performed perturbatively. We obtained results 
in the MS scheme at a renormalization scale fi — 2 GeV, 
and extrapolated these to the continuum limit using sev- 
eral lattice spacings. Our final results for the renormal- 
ized rj? are given in Table II for T = and are shown 
in Fig. 5 for T > 0. Combining the results for and 



The value given in Ref. [24] is \l = 3.15(30) GeV" 2 at fj, = 
1 GeV. We divided this by —1.49(7), running the value to the 
scale fi = 2 GeV and accounting for the different sign convention 
we employ, see the remark after Eq. (6). 



the quark condensates we also determined the magnetic 
susceptibilities Xf, see Eq. (26) for the zero temperature 
values. We found Xf to remain constant within errors as 
the temperature is increased up to T ps 170 MeV. 

We showed furthermore that there is a simple relation 
between the tensor coefficients and the spin contribu- 
tion £ s to the total magnetic susceptibility, see Eq. (4). 
The negative sign of £ s reveals a diamagnetic response, 
i.e., that the spin magnetization of the medium aligns 
itself antiparallel to the external field. The magnitude 
of this effect reduces as the temperature grows, as £ s is 
proportional to rj, plotted in Fig. 5. For the free case £ 
and £ L are known to have opposite signs [52], implying 
a partial cancellation between the two sectors. There- 
fore, a determination of the orbital angular momentum 
contribution is necessary to arrive at a definite conclu- 
sion whether the total response of the QCD vacuum to 
external magnetic fields is para- or diamagnetic. 
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Appendix A: Spin- and orbital angular momentum- 
contributions 



The partition function of QCD (this time without tak- 
ing the root of the determinant as in the staggered lattice 
formulation) is given by the functional integral, 



Z = I VUe~ 0S <> J|det(0 
J f 



TO/), 



(Al) 



with the massless Dirac operator Tp t = j^D^j and co- 
variant derivative D^j = + iqfA^ + iA^T a . For 
an external magnetic field in the z-direction one has 
d x A y - d y A x = B and A z = A t = 0. 

The derivative of the logarithm of Eq. (Al) with re- 
spect to B is 



dlogZ 
dB 



f 



tr- 



1 



dH>, 



dB 



(A2) 



We manipulate this using tr dip t /dB cx tr 7^ = and 
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the cyclicity of the trace: 
d log Z 

' f 



dB 



V — (ii- 



- 1 



f + rrif ) dB 
dlp f 



E m/ \ tr $ f + m f ^ f dB 
1^ 1 /, 1 50' 1 



(A3) 



tr- 



The derivative of the square of the Dirac operator in the 
magnetic field background, after a standard simplifica- 
tion involving 7-matrices, reads 



dip) dB 

~d~B = ~d~B 



f 



Qf a xy, 



(A4) 



where B 2 = D^jD^j with summation for jj, but not for 
/. This implies, 



TdlogZ 1 



V dB 2 m f 



E ((^f^y^f) + (^f L ^f)) 
f 



where we defined 



dB) 

~d(q f B)- 



(A5) 



(A6) 



This operator corresponds to a generalized angular mo- 
mentum, as for the choice A x = —By/2, A y = Bx/2 
(such that d^A^ — 0), it assumes the form L xy = 
-i(xd y - yd x ) + q f B{x 2 + y 2 )/2 - yA%T a + xA a T a . 



Eq. (2), the spin-dependent contribution equals 



d 2 Af & 



d{eBf 



= -N, 



e.B=0 



E 



(q f /e) 2 f d 2 p 



2tt 



1 



(2tt) 2 p 2 



m 



f 

_(A9) 

In Appendix B we will calculate the tensor polarization in 
the free case. Comparing Eq. (A9) with Eq. (B5) below, 
we see that the first term of Eq. (A7) is indeed the spin- 
related contribution, £ s . The second term of Eq. (A7) 
is then identified with the orbital momentum coupling. 
The two contributions to Eq. (3) then read, 



e s = E 



(<If/ef 



2m/ 



qf/ed(i/>fL xy ipf) 



f .... , s ^ 2m. d(eB) ' 

(A10) 

where we used the definition of the tensor coefficient, 
Eq. (5). This shows that the tensor coefficient of the 
quark condensate is responsible for the spin contribution 
of the total magnetic susceptibility. Recalling the rela- 
tion between the sign of £ and para/diamagnetism as 
discussed in the introduction, we conclude that with our 
sign conventions 77 > (xf > 0) corresponds to param- 
agnetism, while Tf < (xj < 0) to diamagnetism. We re- 
mark that on the lattice £ L cannot directly be computed 
from Eq. (A6), due to the quantization of the magnetic 
flux. 



Appendix B: Logarithmic divergence in the tensor 
polarization 



Altogether, using the definition of the (total) magnetic 
susceptibility, Eqs. (2) and (3), we get 

= Qf/e ( d{ip f a xy ip f ) d (4>fL xy ip f ) \ 
* f 2m } \ d{eB) ' d(eB) ) g=Q ' 

(A7) 

showing two separate contributions £ s + £ L to the total 
susceptibility, cf. Eq. (4). 

The conventional calculation of the spin- and orbital 
momentum-related contributions to £ yields the same re- 
sult. Below we demonstrate this for the free case. Here 
the spin-related contribution to the change in the free 
energy density due to the magnetic field at zero temper- 
ature is given by [52, 54, 55], 

^ =~ N >lwf < A8 > 

x (yp 2 +m 2 f + sqfB - ^Jp 2 + m'jj . 

/,s=±l 

Employing the definition of the total susceptibility, 



In this appendix we will demonstrate the appearance 
of a logarithmic divergence in the tensor polarization of 
the condensate. We consider one free quark with electric 
charge qf and mass m/ at vanishing temperature. 

The negative square of the Dirac operator in the back- 
ground of a constant magnetic field is well-known to have 
eigenvalues [14, 56] 

-U>) -^\ 2 =pl+p 2 z + (2n + l)\q f B\ + sq f B, (Bl) 

being twice degenerate (incorporating particle and an- 
tiparticle). Here pq, p z are momenta, n — 0, 1, ... labels 
the Landau levels and s = ±1 is twice the spin (these 
are the eigenvalues of <r xy ). which is coupled to the mag- 
netic field (here we do not consider anomalous magnetic 
moments). The sum over the eigenvalues is performed 
according to (see e.g. Ref. [52]), 

^ = 2Nc 2^f J Po ^ ^'H^E E ■ 

A 2 J ~°° J ~°° n=0 »=±1 

(B2) 

For the tensor polarization of Eq. (10) we note that 
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due to chirality (and since 75 commutes with a xy 



Appendix C: Definition of staggered operators 



1 1 

tr^s a x V = tr7 5 -^ — ; 75 a. 



f + «*/ 



7 f + m f 



■>■!) 



(B3) 



tr- 



1 



a xy = m/tr— -J 



which results in the spectral representation (omitting the 
staggered factor of 1 /4) , 



(ipfPxyipf) = N ( 



m,f\qfB\ f d 2 p 



7T 7 (2tt) 



(B4) 



^ p 2 + (2n + 1 + s sign( g/ B)) 



In the sum the contributions {n = fc, s sign(g/-B) = 1} 
and {n = fc+1, ssign(qyS) = —1} cancel leaving only the 
unpaired lowest Landau level {n = 0, ssign(qfB) = — 1}, 
as was also noted in Ref . [29] . Hence we get 



{ipfOxy^f) = -N c 



rrifqfB f d 2 p 1 



7T J (27r) 2 p 2 + mj 



(B5) 



This cancellation can be confirmed via zeta function reg- 
ularization and is absent for other observables like the 
free energy or the condensate. As the eigenvalue of the 
lowest Landau level is ^-independent, the free tensor po- 
larization is exactly linear in the magnetic field. 



To discretize the continuum tensor bilinear operator on 
the lattice in the staggered formalism, one has to define 
fields that live on corners of four-dimensional elementary 
hypercubes of the lattice. The position of a hypercube 
inside the lattice is denoted by the index y, where y is a 
four- vector with even components y^ = 0,2,..., — 2 
(all lattice extents are even). The position of a fermion 
field component within a specific hypercube is defined by 
one additional four- vector index, C (C„ € {0,1}). Fol- 
lowing the notation of Ref. [44], this hypercube field is 
denoted in this section by x(y)c = x(y+C) /4 (instead of 
ij){x), to distinguish between the different ways of label- 
ing). After rotating into the staggered basis, the operator 
O s = ipfT s ipf can be written as [44] 



CD 



x(y)c (ps 



1 



CD 



x(y) 



D-: 



(CI) 



where the matrix T s acts on spin and the unit matrix 1 
on taste components of the staggered field x(v) ( we wn l 
consider the scalar r§ = 1 and the tensor Tj- = spin 
structures). Here 



1 



7^ r s 7 1 



with the staggered transformation given by 
7c = 7i 7l 72 C2 7 3 C3 74 C4 - 



(C2) 



(C3) 



We evaluate the remaining logarithmically divergent 
integral with dimensional regularization in d = 2 — e di- 



mensions. 



(il) f a xy ip f )=N ( 



4ir 2 




(B6) 

A log m^-term has appeared, whose coefficient is scheme- 
independent, for 3 colors it is 3/(47r 2 ) • rrifqfB (cf. 
Ref. [15] with different sign conventions). Also the singu- 
larity for e —> has been isolated and can be subtracted 
through a particular renormalization scheme, introduc- 
ing a cut-off A such that {^PfUxy^f) oc log(m 2 /A 2 ), or, 
on the lattice \og{m 2 fa 2 ) . The finite term (7 — log(47r) in 
Eq. (B6)) is scheme-dependent (in our lattice scheme it 
reads 0.1549 it 2 — log 4) but, together with the logarithmic 
contribution, it disappears from the combination 



(1 -mfd mf )(i} f a xy tpf) 



m f QfB 
2ir 2 ' 



(B7) 



as we also emphasized in the body of the paper, Eq. (24). 
Note that (1 — rrifd mf ) acting on Eq. (B5) renders the 
integral finite and allows for a direct computation of the 
coefficient of the logarithmic term. 



The operator of Eq. (CI) is clearly not gauge invari- 
ant, since \ an d X are defined in different points of the 
hypercube. To restore gauge invariance, we insert the av- 
erage of products of gauge link variables along all possible 
shortest paths connecting the sites y + C and y + D [44] . 
This average is denoted by Uc.d and the gauge invariant 
operator is now 



Os = E *(»)c (r, ® i) CD u c ,d x(v)d 



(C4) 



C.D 



Using the definition of Eq. (C2) we can further simplify 
the expression for the operator O s . One useful identity 
is 



7m 7c = 7C+A tyi(C) , 

so that 
1. 



(0 = (_i)5W^ (C5) 



Tr 



7c 1 7d 



C.D ■ 



1, 



Tr 



7c °> 7d 



= 6, 



(C6) 



(Here and below, in expressions such as C + p, the sum 
is to be taken modulo 2.) Using Eq. (C6), the scalar and 
tensor operator (the latter for any fi 7^ v) can be written 
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as 

o s = J2x(y)Dx(y)D, (C7) 



Ot = - Y] x(y)c+A+i> C^+a+^-d x(y)c 

1 d (C8) 

x r iv {D)ti VL {D + P). 
As can be seen from the above equation, the tensor op- 
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